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We investigate the effects of exponentially correlated noise on birhythmic van der Pol type oscillators. The 
analytical results are obtained applying the quasi-harmonic assumption to the Langevin equation to derive an 
approximated Fokker-Planck equation. This approach allows to analytically derive the probability distributions 
as well as the activation energies associated to switching between coexisting attractors. The stationary proba¬ 
bility density function of the van der Pol oscillator reveals the influence of the correlation time on the dynamics. 
Stochastic bifurcations are discussed through a qualitative change of the stationary probability distribution, 
which indicates that noise intensity and correlation time can be treated as bifurcation parameters. Comparing 
the analytical and numerical results, we find good agreement both when the frequencies of the attractors are 
about equal or when they are markedly different. 

Keywords .".Stochastic bifurcation; colored noise; birhythmic system; Fokker-Planck equation 
Corresponding author: E-mail address; filatrella@unisannio.it (G. Filatrella). 

To appear in: Commun. Nonlinear Sci. Numer. Simulat. 

PACS numbers: 

05.10.Gg Stochastic analysis methods (Fokker-Planck, Langevin, etc. 

82.40. Bj Oscillations, chaos, and bifurcations 
05.45.-a Nonlinear dynamics and chaos 

82.40. -g Chemical kinetics and reactions: special regimes and techniques 


I. INTRODUCTION 

When studying a phenomenon, one is generally interested 
in overriding effects. Nonlinear science has been used to ac¬ 
count for some phenomena in physics, biology and modeling 
of social facts that could not be otherwise explained, even at 
a qualitative level. Among these phenomena, we can number 
self-sustained oscillations and noise activated transitions. A 
paradigmatic model for oscillating systems is the van der Pol 
circuit, named after the work of Balthazar van der Pol, who 
introduced relaxation oscillations to describe the cycles pro¬ 
duced by self-sustained oscillating systems such as the triode 
circuit.With its many variants, van der Pol oscillators 
have served as a paradigm in various physical, chemical, and 
biological processes A special merit of van der Pol 

has been to propose a dimensionless (or reduced) equation 
in the terminology of Curie for it can be used to ex¬ 

plain various system regardless their origins. The van der Pol 
system was not the first one where self oscillating behaviors 
were observed, but this is the first reduced equation in Curie’s 
sense. With its graphical representation van der Pol also 
contributed to change the way in which these nonlinear phe¬ 
nomena had to be investigated. 

The oscillatory motion can be influenced by internal or ex¬ 
ternal noise In UiH] and the resulting interplay between nonlin¬ 
earities and noise can induce transitions similar to the standard 
transition between two stable points ifl^ [oll . In the case of 
self-sustained oscillations, one observes a transition between 
two stable orbits rather than two points. However, it has been 


shown that the escape is indeed of the Kramers’ type for un¬ 
correlated noise ifldl [fsll even in presence of time delay ifl^ 
in , and it can be ascribed to the existence of a quasipotential 
or pseudopotential ifl^l^ . The quasipotential allows to esti¬ 
mate the energy barriers, and indicates that for most parame¬ 
ters, a single attractor is dominantly stable ifldUfSl] . while the 
other is characterized by a small energy barrier. The effect 
of a finite correlation time has also been investigated, and it 
has been found that, as expected for static potential, also the 
periodic orbits are stabilized by the finite correlation 1^ . 

Noise, however, can produce also effects that are more radi¬ 
cal, for it can induce a structural (or topological) change of the 
solutions ||23] . This is the case when noise induces a stochas¬ 
tic bifurcation, that may be characterized with a qualitative 
change of the stationary probability distribution. Moreover, 
noise can stabilize unstable equilibria and shift bifurcations, 
or it can even induce new stable states that have no determin¬ 
istic counterpart (for instance noise can excite internal modes 
of oscillation and it can even enhance the response of a non¬ 
linear system to external signals ll23l - l^ ). When the noise 
driven stochastic differential equation is colored, the nature 
of the stochastic process becomes non-Markovian and cannot 
be treated with the standard Fokker-Planck techniques. The 
study of stochastic nonlinear systems driven by colored noise 
has been undertaken by a number of investigation ll26l - l^ . 
On a general ground correlated noise appears as the result of 
the gross approximation of a system over some finite scale that 
masks the effect of slow variables 11^ . as it has been argued in 
ecology or electronic circuits i24|] . Birhythmic oscillators 
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are no exception, as thw emerge in biochemical reactions ll^ 
or sleeping patterns Jj, where such coarsening could oc¬ 
cur. In many treatments, one is essentially concerned with an 
equilibrium thermal bath at a finite temperature, which stim¬ 
ulates the reaction coordinate to cross the activation energy 
barrier. In this paper, we wish to investigate the transition be¬ 
tween orbits characterized by different frequencies 

(and hence birhythmic systems). To do so we analytically 
treat the effects of Gaussian correlated noise on a special van 
der Pol oscillator system that displays a birhythmic behavior, 
for the two coexisting attractors characterized by different fre¬ 
quencies. We search the relation between stochastic bifurca¬ 
tions and noise correlation time on the dynamical properties, 
taking the system parameters and the statistical characteristics 
of the noise (e.g., noise intensity and noise correlation time) 
as bifurcation parameters. Our approach of the bifurcations is 
based on the stochastic averaging method jsl, that is separat¬ 
ing fast and slow variables of the oscillator. 


The work is organized as follows. In Sect. |n]we discuss 
the model and the numerical approach to the discretized equa¬ 
tions. In Sect. m we derive an analytic approximation that 
includes the effects of a finite correlation time. Section m 
compares the analytical and numerical results. SectionlVlcon- 
cludes. 


II. THE CORRELATED NOISY MODEL AND 
BIRHYTHMIC PROPERTIES 

In this Section we present the van der Pol model that we 
use as a prototype of birhythmic model. We do so first pre¬ 
senting the electrical equivalent circuit and then deriving the 
normalized equations. 

A. The correlated noisy model 

The model used in our analysis can be reproduced with 
electronic circuits. An example of equivalent circuit is shown 
in Fig. [T] It consists of electronic multipliers Mi (i = 
1,...,4), integrators that are operational amplifiers with a 
feedback capacitor, and sommators realized by operational 
amplifiers with multiple input resistors. Using Millman law, 
the characteristic of each component and the contributions of 
the electrical voltage 14 supplied by an external source (that 
in this study is assumed to be a random term) we find for the 
potential at the left of the point S: 

" Ri2 Rio K Rq Rii K^' ^ 

where AT is a scaling factor, which has the dimension of volt¬ 
age. The multiplicator Mi gives the end voltage (the dots 
denote the time derivative); 
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If we take the time derivative of Eq.(|2|, we get 
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It is convenient to introduce the variables t' = ojot and 
V = Vqx, where Vq is the reference voltage. With the con¬ 
straint Ri 2 Vq = RiqEK, Eq. (|3ll becomes the following 
non-dimensional differential equation; 


with a cubic function truncating the Taylor-Mclaurin series 
expansion Hi. Analogously, one can consider a further ex¬ 
pansion to the 5*^ order 13511^ to retrieve an equation as the 
above @. 


X — /i(l — x"^ + ax'^ — I3x^)x + X = (3) 

where 

R 13 E Ri 2 V^ Ri 2 V§ 

^ KR12R2C1UJ0' “ RgEK^’ ^ RiiEK^’ 


R 2 R 3 CI ’ 




A similar equation was originally obtained by van der Pol in 
the analysis of a triode, when the triode was approximated 


Equation ([3]), even when p(f) 7 ^ 0, is of the van der Pol 
type equation, inasmuch most of the practical applications are 
driven by weak noise intensities. It is also worth mentioning 
that van der Pol provided a stability criterion for periodic so¬ 
lutions when the coefficients are sufficiently small. This equa¬ 
tion is also used to model coherent oscillations in biological 
system H [ 3 ^ [ 3 ^ . In the context of enzymatic reactions lf30l] 
the parameters a and (3 are positive parameters which mea¬ 
sure the degree of tendency of the system to a ferroelectric 
instability compared to its electric resistance, while p is the 
parameter that tunes non linearity. 
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FIG. 1. Circuit diagram for an electronic circuit approximately described by a birhythmic van der Pol type equation. Here the Mi’s are analog 
multipliers circuits while the operation amplifiers with feedback capacitors perform analogue integrations. 


We assume that both the noise r]{t) is stationary and in¬ 
dependent of the memory Kernel, thus one cannot apply the 
fluctuation-dissipation theorem. We also assume that the ran¬ 
dom excitation 77 (f) is a zero average: 

( 77 (f)) = 0 (4) 


The parameter r can be also referred to the intrinsic correla¬ 
tion time: 


r = 


2D 

W) 


( 10 ) 


correlated noise, so we have 

V{i)) = DX exp [-A(f - f)]. (5) 

(the parameter D is the intensity of the noise and r = is 
the correlation time for the colored noise.) 

The correlated noise pit) can be generated as the solution 
of the Langevin equation : 


77 (f) = -Xr]{t) + Xgw{t) ( 6 ) 


where (f) is a standard Gaussian distribution of unit vari¬ 
ance. 

To complete the description, one further specifies the dis¬ 
tribution of initial 770 values in Ea.dTali. denoted by the paren¬ 
thesis {...}:, which is essential for the stationary correlation ( 
one should consider different initial conditions, with different 
initializations of the correlated noise, a caution not necessary 
for uncorrelated noise): 

{(? 7 (f)? 7 (f))} = fAAexp(-A|f-f|). (7a) 

The distribution of the initial values is given by: 


P{Vo) 





( 8 ) 


The overall noise intensity D is the zero-frequency part of the 
power spectrum of the (stationary) noise source 77 (f) 


2D 


/ OO 

|(77(f)y7(0))|df. 

-00 


(9) 



a 


FIG. 2. Parameters domain for the existence of a single limit cycle 
(white area) and three limit cycles (colored area) with p, — 0.01 and 
the parameter Si{a, l3){i = 1, 2, 3, 4, 5, 6) vrAere fli = 123 = 1 
(Table\B and Si{i = 7, 8, 9,10) where !2i / Q 3 (Table\I^. 


The analysis of the noisy dynamical flow involves a study in 
terms of a two-parameter space (D,f). Accurate approxima¬ 
tion schemes for colored noise are only valid in the asymptotic 
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Si{a,p) 

Amplitudes Ai of the orbits 

Frequencies Hi of the orbits 

Periods Pi of the orbits 


Ai=2.3772 

ni=1.0021 

Pi = 6.270 

Si (0.114; 0.003) 

A2=5.0264 

112 = 1.0011 

P2 = 6.275 


As =5.4667 

123=1.0231 

P 3 = 6.141 


Ai =2.3069 

Hi=0.9870 

Pi = 6.366 

S 2 (0.1; 0.002) 

A2=4.8472 

112=1.0001 

P 2 = 6.275 


A3=7.1541 

113=0.97123 

P 3 = 6.468 


Ai =2.4269 

Hi=0.9850 

Pi = 6.379 

S 3 (0.12; 0.003) 

A2=4.2556 

H2=0.9990 

P 2 = 6.289 


A3=6.3245 

113=0.9865 

P 3 = 6.369 


Ai =2.4903 

Hi= 1.0002 

Pi = 6.282 

S4(0.13; 0.004) 

A 2 =4.4721 

H 2 =i.oooi 

P 2 = 6.282 


A3=5.0791 

H3=0.9991 

P 3 = 6.289 


Ai =2.6605 

Hi=i.ooo 2 

Pi = 6.282 

Ss (0.145; 0.005) 

A2=3.8305 

H 2 =i.oooi 

P 2 = 6.282 


A3=4.9643 

113 = 1.0005 

P 3 = 6.280 


Ai=2.7864 

Hi =0.9992 

Pi = 6.288 

S6(0.154; 0.006) 

A2=3.8821 

H2=1.0001 

P 2 = 6.282 


A3=4.2698 

113 = 1.0002 

P 3 = 6.282 


TABLE I. Characteristics of the approximated JU/ limit cycles, as per Eg. jilt, of the birhythmic van der Pol system, EaAllt for the noiseless 
case (D = 0) when the two stable frequencies are about equal ("i.e. ~ QsJ. All data refer to the case p = 0.01. 


Si{a,l3) 

Amplitudes Ai of the orbits 

Frequencies Hi of the orbits 

Periods Pi of the orbits 

S7(0.12; 0.0014) 

Ai= 2.491378 
A2=3.52558 
A3=10.88605 

Hi=1.00210 

H2=0.99994 

H3=0.57300 

Pi = 6.2698 

P 2 = 6.2834 

P 3 = 10.9651 

S8(0.09; 0.001) 

Ai= 2.26969 

A2= 4.59373 

A3= 10.85109 

Hi=0.99986 

Ha =0.99974 

H3=0.68107 

Pi = 6.2839 

Pa = 6.2847 

P 3 = 9.2252 

S9(0.12; 0.0016) 

Ai= 2.48185 

A2= 3.58637 

A3= 10.04878 

Hi=1.00034 

H2=0.99993 

H3= 0.72935 

Pi = 6.2809 

P 2 = 6.2845 

P 3 = 9.2352 

Sio(0.12; 0.0024) 

Ai= 2.44836 

A2= 3.88657 

A3= 7.67464 

Hi=0.98930 

H2=0.99989 

H3=0.86703 

Pi = 6.3510 

P 2 = 6.2837 

P 3 = 7.2466 


TABLE 11. Characteristics of the approximated H^l limit cycles, as per Eq. GH, of the birhythmic van der Pol system, Ea. iTlt for the noiseless 
case (D = 0) when the two stable frequencies are different, ("i.e. Hi 7 ^ CI 3 ). All data refer to the case p = 0.01. 


limits of one the both parameters D and/or r, as discussed in 
Ref. ll^. 

B. Birhythmic properties 

The properties of the dynamical attractors of the modified 
van der Pol model Q can be analyzed as follows. One first 
assumes that the periodic solutions of Eq. @ are represented 

bylS 

x{t) = Ai cos Hit. ( 11 ) 

The amplitudes Ai and the frequencies Hi can be analytically 
approximated JUt]; in this scheme the amplitude Ai is inde¬ 


pendent of the coefficient p, that only enters in the frequency 
Hi of the orbits. The number of cycles of the modified van 
der Pol system Eq. © is determined by the parameters a 
and /3. It is found that in some region of the parameter space 
one obtains three limit cycles (two of them are stable, and 
one is unstable). For each of the two stable limit cycles it 
corresponds to a different frequency, and hence the system 
is birhythmic. The parameter plane (/3-a) is represented in 
FiglU where a shadowed area denotes the portion of the re¬ 
gion where birhythmicity occurs At the border of the 

shadowed area one observes the passage from a single limit 
cycle to three limit cycles. The limit cycles exhibit very sim¬ 
ilar frequencies in the region of low a (Table |I]i and clearly 
different frequencies when a increases (Table HU). By way of 
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conclusion of this part, we have described the equation of mo¬ 
tion of the system of Fig. [T]when the current is truncated to 
the fifth order. The system is known to be birhythmic, and the 
main properties, approximating the solution as in Eg.lfTTIl are 
summarized in TablesHlandHIl 

III. ANALYTIC ESTIMATES 

In this Section we derive the properties of the attractors de¬ 
scribed in Sect. |n]when subject to correlated noise. 

A. Residence times of the attractors 

The evolution of the parameter a; (f) is subjected to the influ¬ 
ence of random forces, which permanently tend to destabilize 
it. Random fluctuations, in the case of three attractors (two 
stable and one unstable) can induce a transition between the 
two stable attractors |fl4|] . The ubiquitous problem of noise 
induced escape from a metastable state is characterized by a 
rate; when the phenomenon occurs on a long time scale com¬ 
pared to dynamic time scales of the localized solution it as¬ 
sumes the features of an Arrhenius like behavior for a dynam¬ 
ical attractor and in particular in van der Pol birhyth¬ 

mic systems lU^fl^ . 

Fig. [3 is an illustration of the escape process from a 
state of local stability via noise-assisted hopping. The stable 
states are characterized by an asymmetric double well pseudo¬ 
potential U (A) with angular frequency of the metastable state 
ujAi 3 and positive-valued angular frequency of the unstable 
state at the barrier ijja 2 defined through the relation ^ = 

and \ [/"(A 2 ) |. Thus, consid¬ 

ering the particle of mass Mg moving in a double well quasi¬ 
potential, the rate coefficients K'^ and K~ are the forward 
and backward rates, respectively: 

= v± exp (12) 

where A(7± denote the threshold energies for activation and 
v± are the corresponding pre-factors. 

The inverse of the Kramers escape rate is the escape time 
from one well to another. The two times might be very differ¬ 
ent 11^1431] . To measure the different properties, we compute 
the average persistence or residence time Ti 3 on the attractor 
with limit cycle amplitude Ai 3 as: 


Ti = T(Ai ^ A 3 ) = 

1/K+ 

(13a) 

1/K+ + \/K-' 

T 3 = r(A3 ^ Ai) = 

1/K- 

(13b) 

1/K+ + 1/K-' 


We wish to emphasize here the special character of both the 
coordinate A and of the quasipotential U , see Fig. [3 In fact A 
is the amplitude of the oscillations approximated by Eq. (HB, 
and therefore it represents the amplitude of an approximated 


orbit, rather that a single point. Also, the quasipotential ?7 is a 

(a) 




FIG. 3. Schematic illustration of the escape process in an asymmetric 
metastable potential. Escapes occur via the forward rate and the 
backward rate K~, respectively, while AUi ,3 are the corresponding 
activation energies. 


Lyapunov function that characterizes the asymptotic behavior 
of the nonequilibrium state ifT^ , not a bona fide potential. 
Therefore Eig. [3 is, at this stage, a pictorial description of 
the process. However, we will show in the next subsection 
that a function playing the role of an activation energy (i.e., 
determining the rate of rare escapes in the low noise intensity 
limit), can be analytically retrieved. 


B. Stationary probability distribution 

In the quasi-harmonic regime, on the assumption that the 
noise intensity is small, it is convenient to use a change of 
variables suggested by the approximation Eg.lfTTIl. i.e. x = 
A cos 9 and 9 = ip + t (we assume H = 1). The instantaneous 
amplitude A{t) and phase (pU) are given by the following ef¬ 
fective Langevin equation 
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a a 

FIG. 4. Ejfect of the noise and of the correlation time r on the boundary between the region of one and three limit-cycle solutions in the 
parametric {a, P)-plane for the stationary probability current governed by Ea. M3^ as per the analytic solution Eq. {Zl. (a) The effect of 
correlation time t for D = 0.1; (b) The effect of noise intensity D for correlation time r = 0.1. 


dA pA 


[64 - 16A^ + 8aA'^ - 5/3^®] + 


D * / X dU ^ , , 

2A(l + r2) = 


here and represent the noise sources, and U is an effective potential or quasipotential; 


U = - 


pA 


32^2 - fA'^ + ^aA^ - ^I3A^ 


D 


2(1+ t2) 


\nA. 


(14a) 

(14b) 


(15) 


Contrary to the case of uncorrelated sources, a rigorous 
derivation of the noise characteristics cannot be obtained, for 
the additional variable rj introduces an extra dimension. How¬ 
ever, dA does not depend on (p, thus we assume that we can 
develop a probability density for A, rather than a joint den¬ 
sity for A and ip. For 7 ^ 0, we also assume that Ea. (ll4bl i 
amounts to an equation for the coordinate A driven by corre¬ 


lated noise and we ignore the fact that it represents an orbit 
rather than a single point. 

Consequently, we treat the effect of noise as in the case 
of a standard double well potential associated to a driven os¬ 
cillator H. In conclusion, the probability density function 
p(A, f|Ao, to) of the averaged amplitude equation is assumed 
to be governed by the condition that the probability flux Sa 
vanishes, i.e.: 


p{A) :Sa = 0^ dMA) - ^ = 0, 


where 


jjA,A ^ D 

" 2(1+ r2) 


D 


Hence, the stationary solution is; 

p{A) = A^Aexp 


( 256^2 - 32A^ -h —aA^ - 5/3A® 
^ 512H V 3 


(16) 


(17a) 

(17b) 


(18) 


where 
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FIG. 5. Stochastic P-Bifurcation in the parameter plane {D, r) in the case where the two stable frequencies are equal, fli = fla. The dotted 
area represents the parameter region where two sable orbits coexist. Parameters of the system are fi = 0.01 and Si{i = 1, 2, 3,4, 5, 6) 


N = 



r l + r 2 

^ 512D 


1 ^ 256^2 


32^4 + 

3 



is the normalization constant. 


(19) 


Equations (118119b show that in the limit r 0 the colored 
noise ri{f) tends to a white noise 12^ . It is important to note 
that the transitions between the unimodal and the bimodal sta¬ 
tionary probability densities are also referi'ed to as the noise 
induced transitions, and stochastic bifurcation discussed here 
is closely connected to the noise induced transition ifistl . The 
probability distribution is in general very asymmetric, there¬ 
fore for most the parameters a or /3 one can localize the prob¬ 
ability function around a single orbit. The peaks of the prob¬ 
ability distribution can be located deriving the logarithm of 
the probability distribution, d [lnp(y4)] /dA = 0. Thus, to de¬ 
termine the extrema of the distribution Eq. (fTSl) amounts to 


determine the roots of the equation: 


. At 513At 




D 


64 


pAlil + r^) 


= 0 


( 20 ) 


Eor D = 0, the amplitude equation (l20l) coincides with the de¬ 
terministic amplitude equation lU^. Generally this equation 
can admit three positives roots lU^fl^ varying the parameters 
a, /3, /i, and D the number of the real roots of Ea. dTOl) change. 
The effects of the parameters in the system Ea. (118119b can be 
seen as a type of a stochastic P-bifurcation 11^ . i.e a sudden 
change of the probability density function. 

In Eig. m we report the influence of the correlation time t 
and the noise intensity D on the birhythmic behavior; in par¬ 
ticular we show in Eig. |4ta) the effect of the correlation time 
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Noise intensity, D 



FIG. 6 . Stochastic P-Bifurcation in the parameter plane [D, r) in the case where the two stable frequencies are different, Hi 7 ^ 0 , 3 . The 
dotted area represents the parameter region where two sable orbits coexist. Parameters of the system are fj, = 0.01 and Si{i = 7, 8 , 9,10) 


for fixed D = 0.1, and in Fig. Efb) the effect of noise inten¬ 
sity for fixed correlation time (r = 0 . 1 ) as predicted by the 
analysis of Ea. (l20l i. From the Figure, it is clear the opposite 
effect of the noise intensity D and of the correlation time r, 
as the region of birhythmicity increases when r increases, and 
decreases when D increases. 

In the same vein to determine the effects of such noise, we 
represent in Figs. 1516 1 the bifurcation diagram for the sys¬ 
tem parameter plane (D, r) respectively in the case where 
the two stable frequencies are equal (fli = and differ¬ 
ent (rii 7 ^ ria). The domain of the existence of three limit 
cycles decreases when the correlated time r and the noise in¬ 
tensity D increase, and disappears when these parameters be¬ 
come large enough. One finds the same phenomena for all the 
set of parameters Si. 


IV. COMPARISON WITH NUMERICAL SIMULATIONS 

To check the validity of the approximations behind the an¬ 
alytic treatment that has led to the probability density func¬ 
tion (fTSl l. we have performed numerical simulations of the 
Langevin dynamics Q with the numerical scheme to be de¬ 
scribed in Subsect. IIV Al 


A. Numerical algorithm for colored noise 

There are several methods and algorithms to solve stochas¬ 
tic differential equations as the implicit midpoint rule 


with Heun and Leapfrog methods or faster numerical algo¬ 
rithms such as the stochastic version of the Runge-Kutta meth¬ 
ods and a quasisymplectic algorithm Issl - liol] . The starting 
point is the Box-Muller algorithm 1511 to generate exponen¬ 
tially correlated colored noise distributed random variable pAt 
from the Gaussian white noise and two random numbers a and 
b which are uniformly distributed on the unit interval [ 0 , 1 ]. 

To simulate the exponentially correlated colored noise, 
Eqs.Q are replaced with the following equations lE^I : 


X = u, ( 21 a) 

it = /r(l — -\- ax'^ — (3x^)u — x + r], ( 21 b) 

77 =-Ap-F Apu,, (21c) 

where is a Gaussian white noise: 

{9w(t)) = 0 

{gn.{t)g^{t')) = 2DS{t-t'). (22) 


Combining Ea. (l21cl l and Eas. (l22l i one gets an exponentially 
correlated colored noise 77 . To numerically solve Eas. (1211 1. we 
use the Box-Mueller algorithm llddTl to generate a Gaussian 
white noise from two random numbers a and b, which are 
uniformly distributed on the unit interval, as in Refs. inlH. 

The equations have been integrated halving the step size un¬ 
til consistent results are obtained (the problem is particularly 
delicate in the proximity of an absorbing barrier (53l)- Fur¬ 
thermore, the procedure has been calibrated with a standard 
activation barrier to retrieve the Kramers escape rate ifl^ , as 
modified by the correlation time r llz^ . In the end, we have 
find that a step size At = 0.0001 is generally appropriated. 
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Amplitude, A 



Amplitude, A 



Amplitude, A 

FIG. 7. Probability distributions of the amplitude A for different val¬ 
ues of the correlation times r when the frequencies of both attractors 
are identical, i.e fli = fis = 1. Solid lines denote the analytical 
results, while dashed lines denote numerical results. Parameters of 
the system refer to Si of Table\I\with p, — 0.01, D — 0.01 



Amplitude, A 



Amplitude, A 



Amplitude, A 



FIG. 8 . Probability distributions of the amplitude A for different val¬ 
ues of the correlation times r when the frequencies of the attractors 
are different, fli 7 ^ fis. Solid lines denote the analytical results, 
while dashed lines denote numerical results. Parameters of the sys¬ 
tem refer to Si ofTable\n\with p = 0.01, D = 0.4 


but in few cases it has been necessary an even smaller step. 
Moreover, we have averaged the results over 200 realizations, 
that ensures converge within few percents. This scheme has 
been employed to check the validity of the approximations 
behind the analytic treatment that has led to the probability 
density function ( fTSl l. 

In Fig. |7] we plot the behavior of the probability distri¬ 
bution p as a function of the amplitude A for low noise in¬ 
tensity D = 0.01 and several values of the correlation time 
T, when the frequencies of both attractors are similar, i.e. 
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n=0.1 



n=0.01 





FIG. 9. Probability distributions versus the amplitude A for different 
values of the coefficient p, when the frequencies of both attractors are 
different, fli 7 ^ fis. Solid lines denote the analytical results, while 
dashed lines denote numerical results. Parameters of the system refer 
to corresponding to Sg of TableURwith D = 0.05 and r = 0.0001 


(i) We begin with a correlation time that is substantially 
zero (t = 0 . 001 ) and Ti is larger than T 3 ; in such con¬ 
ditions the attractor of the limit-cycle amplitude Ai is 
more stable than the limit-cycle amplitude A^. (This 
correspond to the set of parameters 81 , 83 , 35 .) Increas¬ 
ing the correlation time one observes the symmetric 
case Ti ~ T 3 : both attractors are equivalent and the 
system is in a symmetric bistable double well, i.e. it 


~ 173 ~ 1- It appears that the system is more likely 
found at two distinct distances from the origin, the essential 
feature of birhythmicity. In fact, for the parameters here con¬ 
sidered, there are two attractors in the deterministic system, 
see Table U However, as underlined in the previous Section, 
we expect that noise can destroy birhythmicity (see Fig. |4l). 
In fact, for short correlation times, the amplitude has some¬ 
times a single peak (parameters 54,5 , 5 ). We also notice that 
the correlation noise tends to restore the birhythmic behavior, 
and in fact for the largest values of the noise correlation time 
T we always find a two-peak solution. The probability dis¬ 
tribution p is asymmetric for the set of parameters S'i_ 2 , 3 , 4 , 5 - 
Although the correlation time varies from the values r = 0.77 
to r = 1.47, this asymmetric property of the density proba¬ 
bility distribution p persists. In all these configurations, the 
comparison between the analytic and numerical results is ac¬ 
ceptable. 


Fig. [ 8 ] shows the variation of the probability distribution 
density p versus the amplitude A when the frequencies of both 
attractors are different, (i.e. f7i 173, see Table HH i for D = 
0.4 and different values of the correlation times r. Also in this 
case it is evident that the correlation time r tends to stabilize 
a second orbit, and thus gives rise to a birhythmic behavior. 


Finally, in Fig. |9] we show the probability distribution 
function when the frequencies that characterize the attractors 
are different varying the nonlinear parameter p. The Figure 
demonstrates that the analytic and numerical estimates appear 
to be quite close, except at high levels of the parameter p. 


To measure the effect of the correlation time on the stabil¬ 
ity of the two solutions, let us return to Eq. (fOT l that shows 
the connection between the average persistence (or residence) 
times Ti 3 on the attractor with limit cycle amplitude Ai 3 and 
the corresponding escape rates determined by the quasipoten¬ 
tial (fTSl i. The analytical and numerical behaviors of the res¬ 
idence times versus r are reported in Figs. [T 0 ]and[TT] for 
f7i ~ 173 and f7i 7 ^ 173, respectively. In Fig. [TOl three rel¬ 
evant cases can be found for the first occurrence (f7i ~ 173, 
TablelB: 


has approximately the same probability to stay in one 
or the other basins. As the correlation time is further 
increased, one finds the reverse situation: Ti < T 3 , and 
the first attractor becomes less stable. Thus, for larger 
correlation times the system is more likely to stay on 
the limit-cycle attractor A 3 . It is interesting to notice 
that the agreement between the approximated analysis 
and the numerical simulations is satisfying for all val- 
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FIG. 10. Residence times versus the correlation time t for different values of the parameter Si = {cc, P) when the frequency are almost 
identical f fli = fls = Ij. Solid lines denote the analytical results, while dashed lines denote numerical results. Parameters of the system 
refer to Si of Table\I\with p = 0.01, D = 0.01. 


ues of the correlation time, also when the correlation is 
comparable to the period of the cycle, 27r/n ~ 27r. 

(ii) In the second case when the correlation time is in¬ 
creased from zero (r > 0 . 001 ) the attractor orbit 

is first more stable than the attractor orbit Ai (this cor¬ 
respond to the set of parameters S 2 , S 4 ), and becomes 
less stable with the increases of the correlation time t. 
In other words, this is the reverse transition respect to 
the first case (i); 

(iii) In the third case (corresponding to the case Sq), the at¬ 
tractor orbit ^3 is much more stable (T 3 > Tf), as in 
(ii). However, as the correlation time increases, one 
observes the transition between asymmetric and sym¬ 
metric cases, Ti ~ T^, that seems to be the asymptotic 
solution for extremely large correlation times. 

In Figure [TT] we instead investigate the case when the two 
attractors are characterized by different frequencies (Tablellll). 
we find still another case (for the attractors S^j, S%, and S'g ): 
the inner attractor corresponding to the orbit Ax is more stable 
when the correlation time is small, and becomes progressively 
even more stable when the correlation time is increased. 


In conclusion, the effect of the correlation time on the rel¬ 
ative stability of the two attractors seems to be very differen¬ 
tiated across the parameter space. However, despite the very 
different behaviors, the approximated analysis seems to cap¬ 
ture the main physical effects. 

To outline the behavior as a function of the oscillator pa¬ 
rameters, we show in Figs. [T 2 ]the behavior of the activation 
energies for increasing values of the noise correlation time 
(r = 0 . 001 , 1 , TT, 27r) when the parameter a increases form 
0.095 to 0.135 for /3 = 0.002. In this range of parameters, 
the solution exhibits (almost) identical frequencies of the two 
attractors (Hi .3 ~ 1). In this change of the parameter a the 
role of the two activation energies is exchanged; at low a the 
lowest energy is the inner energy, while at the other limit the 
outer activation energy is the lowest. However, the approxi¬ 
mation (flSl l is valid all through the parameter change for cor¬ 
relation times that are as long as the period of the oscillations 
(r ~ 27r/H). 

Figs. [T3]describes the same type of analysis keeping a con¬ 
stant value of the parameter a = 0.12 and varying the pa¬ 
rameter /?. Also in this case the effective trapping energy un¬ 
dergoes a drastic change and the roles of the two barrier are 
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FIG. 11. Residence times versus the correlation time t for different values of the parameter Si = (a; /3) when the frequency are different { 
Hi 7 ^ ris ). Solid lines denote the analytical results, while dashed lines denote numerical results. Parameters of the system refer to Si of Table 
\i^with p = 0.01, D = 0.2. 


1=0.001 T=1 



FIG. 12. Energy barriers AUi ,3 versus a. The lines with symbols refer to the energy barrier of the inner orbit (AUi, left axis), while the 
lines without symbols to the outer orbit activation energy /AUs, right axis). Solid lines denote the estimates as per Ea. MS^ . while dashed lines 
denote numerical results. The other parameters of the system are p — 0.01 and j3 = 0.002. 


exchanged. Thus, the numerical results of Figs. ll2ll.3l confirm V. CONCLUSIONS 

that the analytic approach leading to the quasipotential (fTSl l is 
very accurate in describing the system. 

The main conclusion of this work has been the existence of 
a drastic change (a P-bifurcation ll^ l of the probability dis¬ 
tribution function of a birhythmic van der Pol like oscillator 
subject to exponentially correlated noise. We have particu¬ 
larly investigated the evolution of the modified van der Pol 
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•c=0.001 T=1 



FIG. 13. Energy barriers A?7i,3 versus jS. The lines with symbols refer to the energy barrier of the inner orbit (Ai, left axis), while the lines 
without symbols to the outer orbit activation energy /AU 3 , right axis). Solid lines denote the analytical estimates as per Ea. MS^ . while dashed 
lines denote numerical results. The other parameters of the system are p, = 0.01 and a = 0.12. 


oscillator in the region where birhythmicity occurs. The char¬ 
acteristics of the birhythmic properties are strongly influenced 
by the nonlinear coefficients as well as the noise intensity and 
correlation. We have performed a detailed analysis of the bi¬ 
furcation diagrams in the parameters plane, and it is evident 
that the noise intensity and the correlation time can be treated 
as bifurcation parameters. An approximated solution of the 
Langevin equation based on the Fokker-Planck equation in 
the quasiharmonic regime gives the stationary probabil¬ 
ity density p{A) of the instantaneous amplitude. With this 
approach, a stochastic bifurcation, a qualitative change of the 
stationary solution, occurs when the intensity or correlation of 
the noise is changed. In this system stochastic P-bifurcations 
correspond to the appearance and disappearance of one 
of the maxima of the distribution of the amplitudes p{A). The 
boundary of the existence of multi-limit-cycle solutions, in the 
parametric (a , /3)-plane, is affected in different manners by 
the noise intensity D and the correlation times t; As expected 
for the standard analysis of correlated noise, to increase the 
correlation time amounts to decrease the effect of the noise in¬ 
tensity. The simulations have confirmed that the approximated 
(obtained through stochastic averaging) Fokker-Planck equa¬ 
tion well describes this type of birhythmic oscillators. We 
therefore expect that real systems (either biological or elec¬ 
tronic circuits) governed by birhythmic van der pol type equa¬ 
tion are very sensitive to changes in the correlation time of the 
noise. In fact, the appearance of correlated noise can drasti¬ 
cally change the structure of the attractors, that could possi¬ 
bly be a signature of the presence of significatively correlated 
noise. For instance, the noise correlation in Fig[TT]or Fig. [12] 
entails that the most stable attractors might become another 
one, with a different amplitude (Fig. fTTI i or both different am¬ 
plitude and frequency (Fig|9l parameter set iSio). 

Summing-up, we have found that the quasiharmonic bal¬ 
ance is (surprisingly) effective in predicting the features of 


the full system, as checked with numerical simulations, even 
when the correlation time is comparable to the period of the 
solutions. Fig. [8| (and this is even more surprising) when 
the two solutions are characterized by different frequencies. 
However, we have found (not surprisingly) that when the pa¬ 
rameter p, that tunes the nonlinearity, is increased, the predic¬ 
tions are less accurate, see Fig. |9| It is therefore tempting to 
draw the general conclusion that the quasi-harmonic balance 
is an effective tool for self-oscillatory systems, and even for 
birhythmic systems characterized by two different frequen¬ 
cies, if the nonlinear terms are kept at bay. Of course, no 
matter however relevant it might be considered the van der Pol 
oscillator, one should also bear in mind that the conclusion is 
based on a specific set of simulations on a single system. As 
such, the conclusion should be taken with due care. 
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